A leap-frog discontinuous Galerkin time-domain method of analyzing electromagnetic scattering problems
School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China
† Corresponding author. E-mail:
yangf@uestc.edu.cn
1. IntroductionNowadays, the world has entered into the era of stealth. However, the analysis of the electromagnetic scattering for low-observable targets, which are often multiscale (the electrically large structures having electrically small details), is still a challenging problem for numerical solvers. The typical frequency-domain method, the method of moment (MoM),[1] is a common choice to accurately deal with these problems. However, the MoM method may become computationally inefficient for wideband computations, since each frequency needs a complete resolution of a linear system of equations. The time-domain method is an attractive alternative since it employs a time-marching algorithm that permits one to find the whole frequency-domain behavior with a single simulation.
Among well-known time-domain methods, the finite-difference time-domain (FDTD)[2] method has become very popular for its versatility and power. However, it requires an orthogonal grid (structured mesh), which will need a high discretization density to capture the geometric characteristics of electrically fine structures, yet this can generate a large number of wasted unknowns in the electrically coarse domains. The sub-gridding technique[3] can alleviate this problem, but it will destroy the simple data structure of the standard FDTD scheme and greatly increase computational complexity. To overcome this limitation, the finite element time-domain (FETD) method, which is more flexible in complex geometric modeling by using an unstructured mesh, has been made[4,5] to solve Maxwell’s equations. But this method requires the solution of matrix equations either directly or iteratively at each time step. A multiscale problem usually is discretized into a great number of unknowns, and it will become expensive to perform operations with a very large sparse linear system of equations during time stepping.
The discontinuous Galerkin time-domain (DGTD) methods,[6–8] which are experiencing a fast development, are promising methods of solving multiscale problems since they possess some of the advantages of FDTD and FETD methods. On one hand, the DGTD method has most of the advantages of FDTD, such as spatial explicit algorithm, memory and computational cost only growing linearly with the number of elements, simplicity and easy parallelization.[9] Meanwhile, the perfectly matched layer (PML) truncation technique[10] can also be straightforwardly integrated into DGTD. On the other hand, the DGTD method inherently possesses the capablity of dealing with arbitrarily shaped and inhomogeneously filled objects, which is an important advantage of the conventional FETD method. The main difference is that the solution is allowed to be discontinuous across the boundaries between adjacent elements, which communicate by means of numerical fluxes.[11–13] Owing to this feature, the DGTD method is more flexible. In Ref. [8], a new three-dimensional (3D) continuous–discontinuous Galerkin finite element time-domain is proposed, it possesses the advantages of a reduced number of unknowns for the continuous Galerkin method and the block-diagonal property of the discontinuous Galerkin method. In a word, the DGTD method is much more flexible than the FDTD method in modeling complex structures and much easier than the conventional FETD method in handling multiscale problems.
There are several key steps for constructing a DGTD system: 1) determining which governing equations the DGTD method will be based on; 2) choosing the element shape and corresponding basis functions for the spatial discretization; 3) applying numerical fluxes to interfaces to stitch all elements together; and 4) selecting a time scheme based on properties of a discretized system.[14]
In this paper, we use a globally explicit DGTD method based on the leap-frog (LF) time integration scheme, to calculate the radar cross sections (RCSs) of some targets, especially including the perfect electric conductor (PEC) and dielectric coated NASA almonds. This geometry has been chosen as a challenging example of low-observation target used in the validation of numerical solvers.[15] At the same time, the electromagnetic fields are expanded into the vector basis functions,[16] which can overcome the shortcomings of node basis functions, such as spurious solutions and difficulty in treating conducting and dielectric edges. The uniaxial PML technique[17–20] is applied to the DGTD modeling. Numerical results show that this method can be competitive with MoM in terms of accuracy.
The rest of this paper is organized as follows. The theory and formulations of this method are described in detail in Section 2. It includes the formulations and the leapfrog algorithms in non-PML and PML region of the DGTD method, and also including boundary conditions. In Section 3, the numerical results are presented, which verify the capability and accuracy by three numerical examples. The conclusions are contained in Section 4.
2. Formulation2.1. DG formulation in non-PML regionThere are mainly two kinds of equations for governing transient electromagnetic problems. Mathematically, these governing equations are equivalent, but with different discretization schemes: they differ from each other greatly in numerical properties.[21] One type is the second-order vector wave equation.[16] The reason for choosing this equation is that the first family of Nédélec elements, also known as the edge elements,[11] which are free of spurious mode, can be directly used for spatial discretization. However, it has difficulties in constructing the time-domain PML, which is a popular technique to truncate unbounded regions. Moreover, it has only one variable, unsuitable for building a DGTD system. Because of the implementation of numerical fluxes, a critical step for the DGTD method is based on both 
and 
variables.
To avoid the above difficulties, the other type, which is the first-order Maxwell’s equations with both 
and 
as field variables, is used for the DGTD method.
The DGTD method is based on a finite-element geometrical discretization of the space M into nonoverlapping elements Vm, each of which is bounded by 
where the electromagnetic field and test functions are expanded into a set of vector basis functions (
. Enforcing the residue of Maxwell’s curl equations to be orthogonal to each basis function, we can obtain
Applying the following vector identity and surface divergency theorem yields
Equation (2) can be rewritten as
where
is the unit normal vector located on surface
and points to the outside of
Vm. The terms on the right-hand side (RHS) of Eq. (
5) are integrals over subdomain interfaces. The computation of the RHS for classical FETD is enforcing the tangential component to be continuous at the interface across adjacent element
(the superscript + denotes magnitudes from adjacent elements, and
. However, in the DGTD method, continuous numerical fluxes of the tangential field components
are used instead of
at each side of
. These tangential fields do not coincide with any of the values on any side of
, but depend linearly on them. So there are several different DGTD systems with different numerical fluxes. One commonly used is the Riemann solver,
[11,12] which is a type of upwind numerical flux and is derived from the physical process of wave propagation and reflection between two different media. The upwind numerical flux is adopted in this paper. This numerical flux can introduce artificial dissipation, and this effect can strongly attenuate the spurious modes.
with
where
is the intrinsic wave impedance of the element
m, and
is the intrinsic wave impedance of the adjacent one.
Another popular and even simpler numerical flux is called the central flux[13]
Considering that the spatial part of the fields is expanded within each element in a set of basis functions equal to the set of test functions, using a Feado–Galerkin method
a final semi-discrete algorithm is found as follows:
where the following items should be noted.
(i) em and hm are column vectors varying in time with the discretized electric and magnetic field coefficients in the element m, and 
and 
with the field coefficients of the adjacent elements.
(ii) Jsm is a column vector varying in time of the discretized excitations in the element m
(iii) M is the mass matrix
(iv) K is the stiffness matrix
(v) L are the flux matrices
2.2. Leapfrog DG algorithm in non-PML regionFor a DGTD system, time stepping can be performed element by element rather than solving a huge matrix system as in FETD schemes. This advantage of the DGTD method can save a large number of memories during time stepping, and furthermore, it makes parallel computation straightforward for a DGTD system.
For the time domain integration, several approaches can be chosen. The most commonly employed one is the second-order LF[22] scheme. The basis of the LF scheme is to sample the unknown fields in a staggered way. The first-order time derivative are approximated by central difference. Average approximation is used to estimate the terms with the electric conductivity. For the two extra dissipative terms arising from the upwind flux formulation, we use the backward approximation, since an average would yield a globally implicit scheme. This fact introduces a slightly more restrictive stability condition.[13] Finally, the resulting fully explicit LFDG algorithm becomes
where the expressions for the constants are
2.3. Boundary conditionsThe flux conditions not only serve as connecting adjacent fields but also implementing boundary conditions.
1) PEC boundary conditions on a face of an element can be strictly enforced by setting the tangential electric field to be null, and tangential magnetic field to be continuous.
Perfect magnetic conducting (PMC) conditions are reciprocal of PEC
2) Regarding the absorbing boundary conditions (ABCs), the straightest one is the so-called first-order Silver–Müller (SM-ABC),[22] which can be enforced by setting 
, 
. For the upwind flux, it is directly implemented by setting 
since it is equivalent to the assumption that there is no contribution to the flux from outside the region of solution. For the centered flux, SM-ABC conditions can also be employed.[22] But in this paper, we mainly have implemented the PML ABCs.[19,20] The detailed formulations of the uniaxial perfect matching layer (UPML) will be discussed in Subsection 2.4.
3) Incident wave conditions can also be generated in a straightforward way. Assuming that a known wave is propagating inside a total-field zone (TFZ), while outside it (scattered-field zone (SFZ)) the field is null. If 
and 
denote the wave fields at each point of the TFZ/SFZ interface (as shown in Fig. 1), the flux across the face (lying on the interface) of an element in TFZ needs to be modified
and if the element is in SFZ
2.4. Formulation in PML regionThe PML plays a critical role in transient simulation of unbounded problems because of its ability to absorb the plane waves of all frequencies and incident angles. Several extensions and improvements have been made since it was proposed by Bérenger.[10,18] Here we introduce the UPML, which is regarded as an artificial anisotropic absorbing material, fulfills Maxwell equations, also attenuates the energy inside the PML without reflection. The UPML can be described in the frequency domain by using Maxwell equations
where the expression of the metric tensor is
One can substitute expression (21) into expression (20), and two auxiliary fields 
and 
are introduced. Then, applying the Fourier transform 
, the set of equations for the PML layer for the fields magnitudes 
, 
and the auxiliary fields 
and 
, can be written as
where all vector magnitudes are expressed in the Cartesian basis. The tensors
,
, and
take the following forms:
In this study, 
, where d and 
are the distance and maximum distance to the interface between the PML and non-PML region, respectively. As indicated in Ref. [23], the value of 
is given by
where
is the average element size in the PML region.
2.5. Leapfrog DG algorithm in PML regionBecause the curl terms in Eqs. (23b) and (23d) do not change from regular Maxwell’s equations, the Galerkin procedure can be straightforwardly applied to it. Meanwhile, considering the fact that the auxiliary fields are expanded with the same set of basis functions, and the auxiliary differential equations are tested following the Galerkin procedure, the semi-discrete scheme for the element m located in the PML region can be found:
where
e and
h are the vectors of the discretized electric and magnetic fields;
K and
L are the stiffness matrix and the flux matrices defined in Eqs. (
14) and (
15), respectively;
M is the mass matrix defined in Eq. (
13);
,
, and
, are mass matrices but affected by the tensors defined in Eqs. (
24). Specific forms of the expression are as follows:
The extension of the LF temporal integration scheme to the semi-discrete system of Eq. (25) is straightforward. Making the usual approximations described previously, the fully explicit algorithm for the PML medium can be derived
where
3. Numerical resultsWe have implemented 3D codes with the upwind numerical flux. The study of the stability of the scheme will not be addressed here, and we have derived heuristic estimations[7,22] for the maximum time steps, yielding stable schemes in each case.
In this section, we show some numerical results of scattering problems that will be illustrated to verify the accuracy and capability of the DGTD method. All simulations are performed on a PC computer with an AMD Athlon(tm) II X4 Processor 3.00 GHz and 16 GB of RAM. The incident plane wave is GAUSS pulse and expressed as
with
and
.
In order to verify the validity of the DGTD method, the first example is performed by a simple case of a scattering problem for a PEC sphere, of which the diameter is 
. As shown in Fig. 2, the simulation region is divided into a total-field zone, holding the sphere, and a scattered-field zone. The surface between both regions is used to excite the plane wave by Huygens sources, through the flux terms in a weak way (for details see Section 2). Another surface is used to compute the near-to-far-field transformation and the RCS. Most of the peripheral surface is the ABCs (UPML is used in this paper), which is used to truncate the whole computational domain.
The structures are illuminated with a vertically polarized plane wave (the electric field vector being along the z axis), which propagates along the negative direction of the x axis (
, 
). The calculated results of the monostatic RCS of the PEC sphere at 
, 
from 0 GHz to 3 GHz are shown in Fig. 3. Meanwhile, it also gives the results of the method of moment (MoM), which are obtained from FEKO software, and the step of frequency sweep is 0.1 GHz. Figure 3 shows that the results from DGTD and MoM are in excellent agreement. The copular bistatic RCSs at 1 GHz, 2 GHz, and 3 GHz in the xoy plane (
, 
and xoz plane (
, 
are shown in Figs. 4 and 5, respectively. It is clearly seen that excellent agreement is again found between the DGTD results and the MoM results.
For verifying the capability to handle the dielectric material problems, the RCS of a dielectric sphere with the relative inductivity 
and the diameter 
mm is analyzed. The simulation regions divided and the information about the incident plane wave are the same as those in the first example. The computational results of monostatic RCS from 0 GHz to 3 GHz with DGTD and MoM are shown in Fig. 6, and are in excellent agreement with each other. The copular bistatic RCSs at 1 GHz, 2 GHz, and 3 GHz in the xoy plane and xoz plane are also given in Figs. 7 and 8, respectively and their excellent agreement is also found.
From two cases above, we observe that the validity and capability to handle the PEC and dielectric scattering problems. However, the analysis of the electromagnetic scattering by low-observable targets is a challenging problem for numerical solvers. In the final example, we use the DGTD method to calculate the RCS of a typical low-observable target: the NASA almond. This geometry has been chosen as a challenging example of a low-observable target used in the validation of numerical solvers.
The NASA almond (Fig. 9) is composed of half ellipsoid: 
and 
,
half elliptic ogive:
and
,
where
d=252.37 mm is the length of the structure. It is a complete double-curvature geometry, which has both smoothly and sharply curved zones, as well as a singular point, the ogive vertex. Apart from a PEC case, a coated material case, which is a perfect dielectric with the relative inductivity
and the thickness
d = 5 mm, has been studied.
The structures are also illuminated with a vertically polarized plane wave (the electric field vector being along the z axis) impinging on the almond at the vertex, see Fig. 9 for the details. The simulation results of monostatic RCS with DGTD and MoM for the same frequency range as the above two cases, are given in Fig. 10. Excellent agreement is found both for PEC and coated cases. At the same time, from the figure it also can be seen that the RCS values of the coating dielectric target are larger than those of the PEC target, and increase with the frequency rising. The higher the frequency, the larger the gap of the RCS between the two cases is. When the frequency is 3 GHz, the RCS value of the coating target is about 11 dBsm larger than that of the PEC target.
The copular bistatic RCSs of the PEC almond at 1 GHz, 2 GHz, and 3 GHz in the xoy plane and xoz plane are presented in Figs. 11 and 12, respectively. It is clearly seen that the DGTD results are in good agreement with the MoM results. Moreover, the copular bistatic RCSs of the coated almond between the two methods are also in excellent agreement in the xoy plane and xoz plane at 1 GHz, 2 GHz, and 3 GHz. The details are shown in Figs. 13 and 14.
4. ConclusionsIn this paper, a fully explicit DGTD method is developed to analyze the scattering of some typical targets, especially the NASA almonds. In the DGTD method the tetrahedral elements and their vector basis functions are used to accurately model the computational domain and electromagnetic fields, respectively. In addition, the UPML is employed to truncate the computational domain. A stable time-difference is used to discretize the equations. The RCS of several typical targets are simulated by using the DGTD method. According to the comparisons of the monostatic and bistatic RCS of the PEC and dielectric sphere for the MoM, the validity and capability of the DGTD method are demonstrated. Furthermore, the comparisons of RCS between the DGTD method and MoM for low-observable targets: the NASA almond, are given in the final example. This excellent agreement of results verifies the accuracy and effectiveness of the DGTD method of dealing with the transient multiscale electromagnetic problems.
